home *** CD-ROM | disk | FTP | other *** search
/ Die Speccy' 97 / Die Speccy' 97.iso / amiga_system / the_aminet / util / misc / cookie21.lha / Cookie / r250.c < prev    next >
C/C++ Source or Header  |  1995-08-22  |  8KB  |  207 lines
  1. /******************************************************************************
  2. *  Module:  r250.c   
  3. *  Description: implements R250 random number generator, from S. 
  4. *  Kirkpatrick and E.  Stoll, Journal of Computational Physics, 40, p. 
  5. *  517 (1981).
  6. *  Written by:    W. L. Maier
  7. *
  8. *    METHOD...
  9. *        16 parallel copies of a linear shift register with
  10. *        period 2^250 - 1.  FAR longer period than the usual
  11. *        linear congruent generator, and commonly faster as
  12. *        well.  (For details see the above paper, and the
  13. *        article in DDJ referenced below.)
  14. *
  15. *    HISTORY...
  16. *        Aug 95: Adjusted the module for use with my Amiga version
  17. *            of a certain cookie program. It seems R250 asserts
  18. *            16-bit ints and 32-bit long ints. Changed int to short.
  19. *            Also commented out floating-point version, and added
  20. *            prototype for the seeding function, and made it more
  21. *            ANSI by adding 'void' here and there.
  22. *
  23. *            Only released as part of this cookie program.
  24. *
  25. *            J÷rgen Grahn - Wetterlinsgatan 13e - 521 34 Falk÷ping -
  26. *            Sweden
  27. *
  28. *        Sep 92: Number returned by dr250() is in range [0,1) instead 
  29. *            of [0,1], so for example a random angle in the
  30. *            interval [0, 2*PI) can be calculated
  31. *            conveniently.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  32. *        Aug 92: Initialization is optional.  Default condition is 
  33. *            equivalent to initializing with the seed 12345,
  34. *            so that the sequence of random numbers begins:
  35. *            1173, 53403, 52352, 35341...  (9996 numbers
  36. *            skipped) ...57769, 14511, 46930, 11942, 7978,
  37. *            56163, 46506, 45768, 21162, 43113...  Using ^=
  38. *            operator rather than two separate statements. 
  39. *            Initializing with own linear congruent
  40. *            pseudorandom number generator for portability. 
  41. *            Function prototypes moved to a header file. 
  42. *            Implemented r250n() to generate numbers
  43. *            uniformly distributed in a specific range
  44. *            [0,n), because the common expedient rand()%n is
  45. *            incorrect.  (J. R. Van Zandt <jrv@mbunix.mitre.org>)
  46. *        May 91: Published by W. L. Maier, "A Fast Pseudo Random Number 
  47. *            Generator", Dr. Dobb's Journal #176.
  48. ******************************************************************************/
  49.  
  50. #include <stdlib.h>
  51. #include "r250.h"
  52.  
  53. /**** Static variables ****/
  54. static short r250_index = 0;
  55. static unsigned short r250_buffer[250] = {
  56.     15301,57764,10921,56345,19316,43154,54727,49252,32360,49582,
  57.     26124,25833,34404,11030,26232,13965,16051,63635,55860,5184,
  58.     15931,39782,16845,11371,38624,10328,9139,1684,48668,59388,
  59.     13297,1364,56028,15687,63279,27771,5277,44628,31973,46977,
  60.     16327,23408,36065,52272,33610,61549,58364,3472,21367,56357,
  61.     56345,54035,7712,55884,39774,10241,50164,47995,1718,46887,
  62.     47892,6010,29575,54972,30458,21966,54449,10387,4492,644,
  63.     57031,41607,61820,54588,40849,54052,59875,43128,50370,44691,
  64.     286,12071,3574,61384,15592,45677,9711,23022,35256,45493,
  65.     48913,146,9053,5881,36635,43280,53464,8529,34344,64955,
  66.     38266,12730,101,16208,12607,58921,22036,8221,31337,11984,
  67.     20290,26734,19552,48,31940,43448,34762,53344,60664,12809,
  68.     57318,17436,44730,19375,30,17425,14117,5416,23853,55783,
  69.     57995,32074,26526,2192,11447,11,53446,35152,64610,64883,
  70.     26899,25357,7667,3577,39414,51161,4,58427,57342,58557,
  71.     53233,1066,29237,36808,19370,17493,37568,3,61468,38876,
  72.     17586,64937,21716,56472,58160,44955,55221,63880,1,32200,
  73.     62066,22911,24090,10438,40783,36364,14999,2489,43284,9898,
  74.     39612,9245,593,34857,41054,30162,65497,53340,27209,45417,
  75.     37497,4612,58397,52910,56313,62716,22377,40310,15190,34471,
  76.     64005,18090,11326,50839,62901,59284,5580,15231,9467,13161,
  77.     58500,7259,317,50968,2962,23006,32280,6994,18751,5148,
  78.     52739,49370,51892,18552,52264,54031,2804,17360,1919,19639,
  79.     2323,9448,43821,11022,45500,31509,49180,35598,38883,19754,
  80.     987,11521,55494,38056,20664,2629,50986,31009,54043,59743
  81.     };
  82.  
  83. static unsigned short myrand(void);
  84. static void mysrand(unsigned short newseed);
  85.  
  86. /**** Function: r250_init  
  87.     Description: initializes r250 random number generator. ****/
  88. void r250_init(short seed)
  89. {
  90. /*---------------------------------------------------------------------------*/
  91.     short          j, k;
  92.     unsigned short mask;
  93.     unsigned short msb;
  94. /*---------------------------------------------------------------------------*/
  95.     mysrand(seed);
  96.     r250_index = 0;
  97.     for (j = 0; j < 250; j++)     /* Fill the r250 buffer with 15-bit values */
  98.         r250_buffer[j] = myrand();
  99.     for (j = 0; j < 250; j++)     /* Set some of the MS bits to 1 */
  100.         if (myrand() > 16384)
  101.             r250_buffer[j] |= 0x8000;
  102.     msb = 0x8000;       /* To turn on the diagonal bit   */
  103.     mask = 0xffff;      /* To turn off the leftmost bits */
  104.     for (j = 0; j < 16; j++)
  105.         {
  106.         k = 11 * j + 3;             /* Select a word to operate on        */
  107.         r250_buffer[k] &= mask;     /* Turn off bits left of the diagonal */
  108.         r250_buffer[k] |= msb;      /* Turn on the diagonal bit           */
  109.         mask >>= 1;
  110.         msb >>= 1;
  111.         }
  112. }
  113.  
  114. /**** Function: r250 Description: returns a random unsigned integer k
  115.                 uniformly distributed in the interval 0 <= k < 65536.  ****/
  116. unsigned short r250(void)
  117. {
  118. /*---------------------------------------------------------------------------*/
  119.     register short    j;
  120.     register unsigned short new_rand;
  121. /*---------------------------------------------------------------------------*/
  122.     if (r250_index >= 147)
  123.         j = r250_index - 147;      /* Wrap pointer around */
  124.     else
  125.         j = r250_index + 103;
  126.  
  127.     new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  128.  
  129.     if (r250_index >= 249)      /* Increment pointer for next time */
  130.         r250_index = 0;
  131.     else
  132.         r250_index++;
  133.  
  134.     return new_rand;
  135. }
  136.  
  137. /**** Function: r250n Description: returns a random unsigned integer k
  138.                     uniformly distributed in the interval 0 <= k < n ****/
  139. unsigned short r250n(unsigned short n)
  140. {
  141. /*---------------------------------------------------------------------------*/
  142.     register short    j;
  143.     register unsigned short new_rand, limit;
  144. /*---------------------------------------------------------------------------*/
  145.     limit = (65535U/n)*n;
  146.     do 
  147.         {new_rand = r250();
  148.         if (r250_index >= 147)
  149.             j = r250_index - 147;      /* Wrap pointer around */
  150.         else
  151.             j = r250_index + 103;
  152.  
  153.         new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  154.  
  155.         if (r250_index >= 249)      /* Increment pointer for next time */
  156.             r250_index = 0;
  157.         else
  158.             r250_index++;
  159.         } while(new_rand >= limit);
  160.     return new_rand%n;
  161. }
  162.  
  163. #if 0
  164. /**** Function: dr250 
  165.         Description: returns a random double z in range 0 <= z < 1.  ****/
  166. double dr250()
  167. {
  168. /*---------------------------------------------------------------------------*/
  169.     register int    j;
  170.     register unsigned int new_rand;
  171. /*---------------------------------------------------------------------------*/
  172.     if (r250_index >= 147)
  173.         j = r250_index - 147;     /* Wrap pointer around */
  174.     else
  175.         j = r250_index + 103;
  176.  
  177.     new_rand = r250_buffer[r250_index] ^= r250_buffer[j];
  178.  
  179.     if (r250_index >= 249)      /* Increment pointer for next time */
  180.         r250_index = 0;
  181.     else
  182.         r250_index++;
  183.     return new_rand / 65536.;   /* Return a number in [0.0 to 1.0) */
  184. }
  185. #endif
  186.  
  187. /***  linear congruent pseudorandom number generator for initialization ***/
  188.  
  189. static unsigned long seed=1;
  190.  
  191.     /*    Return a pseudorandom number in the interval 0 <= n < 32768.
  192.         This produces the following sequence of pseudorandom 
  193.         numbers:
  194.         346, 130, 10982, 1090...  (9996 numbers skipped) ...23369,
  195.         2020, 5703, 12762, 10828, 16252, 28648, 27041, 23444, 6604... 
  196.     */ 
  197. static unsigned short myrand(void)
  198. {
  199.     seed = seed*0x015a4e35L + 1;
  200.     return (seed>>16)&0x7fff;
  201. }
  202.  
  203.     /*    Initialize above linear congruent pseudorandom number generator */
  204. static void mysrand(unsigned short newseed)
  205. {    seed = newseed;
  206. }
  207.